home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / erfc.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  11.3 KB  |  432 lines
  1. /* specfunc/erfc.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  J. Theiler (modifications by G. Jungman) */
  21.  
  22. /*
  23.  * See Hart et al, Computer Approximations, John Wiley and Sons, New York (1968)
  24.  * (This applies only to the erfc8 stuff, which is the part
  25.  *  of the original code that survives. I have replaced much of
  26.  *  the other stuff with Chebyshev fits. These are simpler and
  27.  *  more precise than the original approximations. [GJ])
  28.  */
  29. #include <config.h>
  30. #include <gsl/gsl_math.h>
  31. #include <gsl/gsl_errno.h>
  32. #include <gsl/gsl_sf_erf.h>
  33.  
  34. #include "check.h"
  35.  
  36. #include "chebyshev.h"
  37. #include "cheb_eval.c"
  38.  
  39. #define LogRootPi_  0.57236494292470008706
  40.  
  41.  
  42. static double erfc8_sum(double x)
  43. {
  44.   /* estimates erfc(x) valid for 8 < x < 100 */
  45.   /* This is based on index 5725 in Hart et al */
  46.  
  47.   static double P[] = {
  48.       2.97886562639399288862,
  49.       7.409740605964741794425,
  50.       6.1602098531096305440906,
  51.       5.019049726784267463450058,
  52.       1.275366644729965952479585264,
  53.       0.5641895835477550741253201704
  54.   };
  55.   static double Q[] = {
  56.       3.3690752069827527677,
  57.       9.608965327192787870698,
  58.       17.08144074746600431571095,
  59.       12.0489519278551290360340491,
  60.       9.396034016235054150430579648,
  61.       2.260528520767326969591866945,
  62.       1.0
  63.   };
  64.   double num=0.0, den=0.0;
  65.   int i;
  66.  
  67.   num = P[5];
  68.   for (i=4; i>=0; --i) {
  69.       num = x*num + P[i];
  70.   }
  71.   den = Q[6];
  72.   for (i=5; i>=0; --i) {
  73.       den = x*den + Q[i];
  74.   }
  75.  
  76.   return num/den;
  77. }
  78.  
  79. inline
  80. static double erfc8(double x)
  81. {
  82.   double e;
  83.   e = erfc8_sum(x);
  84.   e *= exp(-x*x);
  85.   return e;
  86. }
  87.  
  88. inline
  89. static double log_erfc8(double x)
  90. {
  91.   double e;
  92.   e = erfc8_sum(x);
  93.   e = log(e) - x*x;
  94.   return e;
  95. }
  96.  
  97. #if 0
  98. /* Abramowitz+Stegun, 7.2.14 */
  99. static double erfcasympsum(double x)
  100. {
  101.   int i;
  102.   double e = 1.;
  103.   double coef = 1.;
  104.   for (i=1; i<5; ++i) {
  105.     /* coef *= -(2*i-1)/(2*x*x); ??? [GJ] */
  106.     coef *= -(2*i+1)/(i*(4*x*x*x*x));
  107.     e += coef;
  108.     /*
  109.     if (fabs(coef) < 1.0e-15) break;
  110.     if (fabs(coef) > 1.0e10) break;
  111.     
  112.     [GJ]: These tests are not useful. This function is only
  113.     used below. Took them out; they gum up the pipeline.
  114.     */
  115.   }
  116.   return e;
  117. }
  118. #endif /* 0 */
  119.  
  120.  
  121. /* Abramowitz+Stegun, 7.1.5 */
  122. static int erfseries(double x, gsl_sf_result * result)
  123. {
  124.   double coef = x;
  125.   double e    = coef;
  126.   double del;
  127.   int k;
  128.   for (k=1; k<30; ++k) {
  129.     coef *= -x*x/k;
  130.     del   = coef/(2.0*k+1.0);
  131.     e += del;
  132.   }
  133.   result->val = 2.0 / M_SQRTPI * e;
  134.   result->err = 2.0 / M_SQRTPI * (fabs(del) + GSL_DBL_EPSILON);
  135.   return GSL_SUCCESS;
  136. }
  137.  
  138.  
  139. /* Chebyshev fit for erfc((t+1)/2), -1 < t < 1
  140.  */
  141. static double erfc_xlt1_data[20] = {
  142.   1.06073416421769980345174155056,
  143.  -0.42582445804381043569204735291,
  144.   0.04955262679620434040357683080,
  145.   0.00449293488768382749558001242,
  146.  -0.00129194104658496953494224761,
  147.  -0.00001836389292149396270416979,
  148.   0.00002211114704099526291538556,
  149.  -5.23337485234257134673693179020e-7,
  150.  -2.78184788833537885382530989578e-7,
  151.   1.41158092748813114560316684249e-8,
  152.   2.72571296330561699984539141865e-9,
  153.  -2.06343904872070629406401492476e-10,
  154.  -2.14273991996785367924201401812e-11,
  155.   2.22990255539358204580285098119e-12,
  156.   1.36250074650698280575807934155e-13,
  157.  -1.95144010922293091898995913038e-14,
  158.  -6.85627169231704599442806370690e-16,
  159.   1.44506492869699938239521607493e-16,
  160.   2.45935306460536488037576200030e-18,
  161.  -9.29599561220523396007359328540e-19
  162. };
  163. static cheb_series erfc_xlt1_cs = {
  164.   erfc_xlt1_data,
  165.   19,
  166.   -1, 1,
  167.   12
  168. };
  169.  
  170. /* Chebyshev fit for erfc(x) exp(x^2), 1 < x < 5, x = 2t + 3, -1 < t < 1
  171.  */
  172. static double erfc_x15_data[25] = {
  173.   0.44045832024338111077637466616,
  174.  -0.143958836762168335790826895326,
  175.   0.044786499817939267247056666937,
  176.  -0.013343124200271211203618353102,
  177.   0.003824682739750469767692372556,
  178.  -0.001058699227195126547306482530,
  179.   0.000283859419210073742736310108,
  180.  -0.000073906170662206760483959432,
  181.   0.000018725312521489179015872934,
  182.  -4.62530981164919445131297264430e-6,
  183.   1.11558657244432857487884006422e-6,
  184.  -2.63098662650834130067808832725e-7,
  185.   6.07462122724551777372119408710e-8,
  186.  -1.37460865539865444777251011793e-8,
  187.   3.05157051905475145520096717210e-9,
  188.  -6.65174789720310713757307724790e-10,
  189.   1.42483346273207784489792999706e-10,
  190.  -3.00141127395323902092018744545e-11,
  191.   6.22171792645348091472914001250e-12,
  192.  -1.26994639225668496876152836555e-12,
  193.   2.55385883033257575402681845385e-13,
  194.  -5.06258237507038698392265499770e-14,
  195.   9.89705409478327321641264227110e-15,
  196.  -1.90685978789192181051961024995e-15,
  197.   3.50826648032737849245113757340e-16
  198. };
  199. static cheb_series erfc_x15_cs = {
  200.   erfc_x15_data,
  201.   24,
  202.   -1, 1,
  203.   16
  204. };
  205.  
  206. /* Chebyshev fit for erfc(x) x exp(x^2), 5 < x < 10, x = (5t + 15)/2, -1 < t < 1
  207.  */
  208. static double erfc_x510_data[20] = {
  209.   1.11684990123545698684297865808,
  210.   0.003736240359381998520654927536,
  211.  -0.000916623948045470238763619870,
  212.   0.000199094325044940833965078819,
  213.  -0.000040276384918650072591781859,
  214.   7.76515264697061049477127605790e-6,
  215.  -1.44464794206689070402099225301e-6,
  216.   2.61311930343463958393485241947e-7,
  217.  -4.61833026634844152345304095560e-8,
  218.   8.00253111512943601598732144340e-9,
  219.  -1.36291114862793031395712122089e-9,
  220.   2.28570483090160869607683087722e-10,
  221.  -3.78022521563251805044056974560e-11,
  222.   6.17253683874528285729910462130e-12,
  223.  -9.96019290955316888445830597430e-13,
  224.   1.58953143706980770269506726000e-13,
  225.  -2.51045971047162509999527428316e-14,
  226.   3.92607828989125810013581287560e-15,
  227.  -6.07970619384160374392535453420e-16,
  228.   9.12600607264794717315507477670e-17
  229. };
  230. static cheb_series erfc_x510_cs = {
  231.   erfc_x510_data,
  232.   19,
  233.   -1, 1,
  234.   12
  235. };
  236.  
  237. #if 0
  238. inline
  239. static double
  240. erfc_asymptotic(double x)
  241. {
  242.   return exp(-x*x)/x * erfcasympsum(x) / M_SQRTPI;
  243. }
  244. inline
  245. static double
  246. log_erfc_asymptotic(double x)
  247. {
  248.   return log(erfcasympsum(x)/x) - x*x - LogRootPi_;
  249. }
  250. #endif /* 0 */
  251.  
  252.  
  253. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  254.  
  255. int gsl_sf_erfc_e(double x, gsl_sf_result * result)
  256. {
  257.   const double ax = fabs(x);
  258.   double e_val, e_err;
  259.  
  260.   /* CHECK_POINTER(result) */
  261.  
  262.   if(ax <= 1.0) {
  263.     double t = 2.0*ax - 1.0;
  264.     gsl_sf_result c;
  265.     cheb_eval_e(&erfc_xlt1_cs, t, &c);
  266.     e_val = c.val;
  267.     e_err = c.err;
  268.   }
  269.   else if(ax <= 5.0) {
  270.     double ex2 = exp(-x*x);
  271.     double t = 0.5*(ax-3.0);
  272.     gsl_sf_result c;
  273.     cheb_eval_e(&erfc_x15_cs, t, &c);
  274.     e_val = ex2 * c.val;
  275.     e_err = ex2 * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON);
  276.   }
  277.   else if(ax < 10.0) {
  278.     double exterm = exp(-x*x) / ax;
  279.     double t = (2.0*ax - 15.0)/5.0;
  280.     gsl_sf_result c;
  281.     cheb_eval_e(&erfc_x510_cs, t, &c);
  282.     e_val = exterm * c.val;
  283.     e_err = exterm * (c.err + 2.0*fabs(x)*GSL_DBL_EPSILON + GSL_DBL_EPSILON);
  284.   }
  285.   else {
  286.     e_val = erfc8(ax);
  287.     e_err = (x*x + 1.0) * GSL_DBL_EPSILON * fabs(e_val);
  288.   }
  289.  
  290.   if(x < 0.0) {
  291.     result->val  = 2.0 - e_val;
  292.     result->err  = e_err;
  293.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  294.   }
  295.   else {
  296.     result->val  = e_val;
  297.     result->err  = e_err;
  298.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  299.   }
  300.  
  301.   return GSL_SUCCESS;
  302. }
  303.  
  304.  
  305. int gsl_sf_log_erfc_e(double x, gsl_sf_result * result)
  306. {
  307.   /* CHECK_POINTER(result) */
  308.  
  309.   if(x*x < 10.0*GSL_ROOT6_DBL_EPSILON) {
  310.     const double y = x / M_SQRTPI;
  311.     /* series for -1/2 Log[Erfc[Sqrt[Pi] y]] */
  312.     const double c3 = (4.0 - M_PI)/3.0;
  313.     const double c4 = 2.0*(1.0 - M_PI/3.0);
  314.     const double c5 = -0.001829764677455021;  /* (96.0 - 40.0*M_PI + 3.0*M_PI*M_PI)/30.0  */
  315.     const double c6 =  0.02629651521057465;   /* 2.0*(120.0 - 60.0*M_PI + 7.0*M_PI*M_PI)/45.0 */
  316.     const double c7 = -0.01621575378835404;
  317.     const double c8 =  0.00125993961762116;
  318.     const double c9 =  0.00556964649138;
  319.     const double c10 = -0.0045563339802;
  320.     const double c11 =  0.0009461589032;
  321.     const double c12 =  0.0013200243174;
  322.     const double c13 = -0.00142906;
  323.     const double c14 =  0.00048204;
  324.     double series = c8 + y*(c9 + y*(c10 + y*(c11 + y*(c12 + y*(c13 + c14*y)))));
  325.     series = y*(1.0 + y*(1.0 + y*(c3 + y*(c4 + y*(c5 + y*(c6 + y*(c7 + y*series)))))));
  326.     result->val = -2.0 * series;
  327.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  328.     return GSL_SUCCESS;
  329.   }
  330.   /*
  331.   don't like use of log1p(); added above series stuff for small x instead, should be ok [GJ]
  332.   else if (fabs(x) < 1.0) {
  333.     gsl_sf_result result_erf;
  334.     gsl_sf_erf_e(x, &result_erf);
  335.     result->val  = log1p(-result_erf.val);
  336.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  337.     return GSL_SUCCESS;
  338.   }
  339.   */
  340.   else if(x > 8.0) {
  341.     result->val = log_erfc8(x);
  342.     result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  343.     return GSL_SUCCESS;
  344.   }
  345.   else {
  346.     gsl_sf_result result_erfc;
  347.     gsl_sf_erfc_e(x, &result_erfc);
  348.     result->val  = log(result_erfc.val);
  349.     result->err  = fabs(result_erfc.err / result_erfc.val);
  350.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  351.     return GSL_SUCCESS;
  352.   }
  353. }
  354.  
  355.  
  356. int gsl_sf_erf_e(double x, gsl_sf_result * result)
  357. {
  358.   /* CHECK_POINTER(result) */
  359.  
  360.   if(fabs(x) < 1.0) {
  361.     return erfseries(x, result);
  362.   }
  363.   else {
  364.     gsl_sf_result result_erfc;
  365.     gsl_sf_erfc_e(x, &result_erfc);
  366.     result->val  = 1.0 - result_erfc.val;
  367.     result->err  = result_erfc.err;
  368.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  369.     return GSL_SUCCESS;
  370.   }
  371. }
  372.  
  373.  
  374. int gsl_sf_erf_Z_e(double x, gsl_sf_result * result)
  375. {
  376.   /* CHECK_POINTER(result) */
  377.  
  378.   {
  379.     const double ex2 = exp(-x*x/2.0);
  380.     result->val  = ex2 / (M_SQRT2 * M_SQRTPI);
  381.     result->err  = fabs(x * result->val) * GSL_DBL_EPSILON;
  382.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  383.     CHECK_UNDERFLOW(result);
  384.     return GSL_SUCCESS;
  385.   }
  386. }
  387.  
  388.  
  389. int gsl_sf_erf_Q_e(double x, gsl_sf_result * result)
  390. {
  391.   /* CHECK_POINTER(result) */
  392.  
  393.   {
  394.     gsl_sf_result result_erfc;
  395.     int stat = gsl_sf_erfc_e(x/M_SQRT2, &result_erfc);
  396.     result->val  = 0.5 * result_erfc.val;
  397.     result->err  = 0.5 * result_erfc.err;
  398.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  399.     return stat;
  400.   }
  401. }
  402.  
  403.  
  404. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  405.  
  406. #include "eval.h"
  407.  
  408. double gsl_sf_erfc(double x)
  409. {
  410.   EVAL_RESULT(gsl_sf_erfc_e(x, &result));
  411. }
  412.  
  413. double gsl_sf_log_erfc(double x)
  414. {
  415.   EVAL_RESULT(gsl_sf_log_erfc_e(x, &result));
  416. }
  417.  
  418. double gsl_sf_erf(double x)
  419. {
  420.   EVAL_RESULT(gsl_sf_erf_e(x, &result));
  421. }
  422.  
  423. double gsl_sf_erf_Z(double x)
  424. {
  425.   EVAL_RESULT(gsl_sf_erf_Z_e(x, &result));
  426. }
  427.  
  428. double gsl_sf_erf_Q(double x)
  429. {
  430.   EVAL_RESULT(gsl_sf_erf_Q_e(x, &result));
  431. }
  432.